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We study earthquake interval time statistics, paying special attention to inter-occurrence times in 
the two-dimensional (2D) stick-slip (block-slider) model. Inter-occurrence times are the time interval 
between successive earthquakes on all faults in a region. We select stiffness and friction parameters as 
tunable parameters because these physical quantities are considered as essential factors in describing 
fault dynamics. It is found that inter-occurrence time statistics depend on the parameters. Varying 
stiffness and friction parameters systematically, we optimize these parameters so as to reproduce the 
inter-occurrence time statistics in natural seismicity. For an optimal case, earthquakes produced by 
the model obey the Gutenberg-Richter law, which states that the magnitude- frequency distribution 
exhibits the power law with an exponent approximately unity. 

PACS numbers: 



I. INTRODUCTION 

Earthquakes are caused by a fracture and frictional slip process. We can understand qualitatively how an earthquake 
occurs on the basis of the plate tectonic theory proposed by A. Wegener. Statistical properties of earthquakes are 
well known as empirical laws while the source mechanism of earthquakes is still an open question. For example, 
Gutenberg and Richtcr proposed the relation between magnitude (M) and frequency (n) expressed by 

logio = a - (1) 

where a and b are positive constants. This relation is called the Gutenberg-Richter (GR) lawj2|. b is the so-called 
&-value and similar to unity. Strictly speaking, b depends on fault structures and seismicity and ranges from 0.8 to 
1.2 fl. 

In general, earthquakes can be categorized into three types: foreshocks, mainshocks, and aftershocks. A mainshock is 
a large earthquake, whereas a foreshock (aftershock) is an earthquake before (after) the mainshock and which occurred 
near the mainshock epicenter. Aftershocks obey the Omori law J|, which stresses that the decay rate of aftershocks 
follows the power law. Subsequently, a modified version was proposed by Utsu Q]. Since the Gutenberg-Richter law 
and the Omori law exhibit the power law, earthquakes are seemed to be self-organized critical phenomena Q|. 

The time intervals between earthquakes can be classified into two types: recurrence times and inter-occurrence times. 
Recurrence time is the interval of time between earthquakes on a single fault or segment, whereas inter-occurrence 
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FIG. 1: Seismicity map around the Japan from 2001/01/01 to 2006/08/31 for M > 5.0. This map is programed on the basis 
of the JMA earthquake catalogs. 



time is the time interval between earthquakes on all faults in a region. For inter-occurrence time statistics, probability 
distributions have been studied by different authors 3, [l^ by using earthquake catalogs (see fig. H]). Recurrence 
times arc generally used by seismologists to describe the time interval between characteristic earthquakes A 



characteristic earthquake is a large earthquake happening on a single fault and depending on fault length, crust 

proposed, such as the log 



12| . However, we cannot 



structure, and so on. For recurrence time, several probability distributions have been 
normal distribution, the WeibuU distribution, and the double exponential distribution 
decide which distribution is the best owing to the lack of data. 

In this work, we focus on the inter-occurrence time statistics. Statistical properties based on earthquake models have 
been investigated numerically and compared with seismicity in nature. Then, earthquake models have been modified so 
as to reproduce fault systems . Generally, numerical simulations have the advantage of enough earthquakes having 
occurred to guarantee statistical accuracy. In addition, it is easy to study the statistical properties under various 
crust conditions by changing control parameters. Optimizing or restricting control parameters so as to adequately 
reproduce the statistical properties, we suggest the probability distribution function of the recurrence time and offer 
new insights into earthquake statistics. 



The stick-slip model proposed by Burridgc and KnopofF [l^ is often called the block-slider model or the Burridge- 
Knopoff model. This model describes the relative motion of faults. Although the model is highly simplified, it has been 
shown that it can extract the statistical properties of earthquakes, such as the GR law, the Omori law, the empirical 



law of the stress distribution, the constant stress drop, and the inter-occurrence time statistics 
This model has been modified in order to describe real crust structures 
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statistics of time intervals in the stick-slip model has not been discussed fully. The purpose of this work is to reveal 
whether the 2D stick-slip model can be understood as useful model in view of the inter-occurrence time statistics. 
Thus, we report numerical investigation of earthquake inter-occurrence time statistics produced by the 2D stick-slip 
model. In this work, stiffness and friction parameters are set as control parameters. Then we restrict or optimize 
these parameters so as to reproduce the inter-occurrence time statistics in nature. It is concluded that the model 
reproduces the inter-occurrence time statistics as well as the GR law in a limited parameter regime. Investigating the 
2D stick-slip model in a optimal case further, we may propose new findings concerning the statistical properties of 
earthquakes. 



II. TWO-DIMENSIONAL STICK-SLIP MODEL 

In this work, we numerically investigate statistical properties of the inter-occurrence time, produced by the two- 
dimensional (2D) stick-slip model. As shown in Fig. 2 (a), the model is composed of blocks on a square lattice, two 
plates, and two different kinds of springs. The upper plate is fixed, whereas the bottom plate moves at a velocity of v. 
So, this model represents the relative motion of faults. A block corresponds to a segment of a fault, so that we define 
a one-block-slip event as a minimum earthquake. Shear stresses and compression stress are modeled respectively by 
coil springs, k^,kj! interconnected by a block, and by the leaf springs, kp between a block and the fixed plate. We 
set the —J/ axis as the direction of the loading plate, and the x axis as the perpendicular to the y axis. Assuming 
that the slip direction of a block is restricted by the y-direction only, the block is described by a stick-slip motion. 
The stick-slip motion can be divided in two parts: one is a stick-state and the other is a slip-state. In the case of the 
stick-state or loading-state, all blocks and the loading plate move together, whereas a block slips on the bottom plate. 
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FIG. 2: (a) 2D stick-slip model, k^, k^, and kp are spring constants. The friction force acts on the surface between the 
block and the bottom plate, (b) a-dependence of the non-linear dynamical friction function, a is fixed at 0.01 throughout the 
simulation. 
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the slip state. 

The equation of motion in a scaled form at cite can be expressed by 



where f/, t', and <f> correspond respectively to a normalized displacement, time, and a dynamical friction force which 
is a function of block velocity. Additionally, 

where v is a maximum of the slipping velocity and vi is the characteristic velocity. The blocks are subject to the 
friction acting on the surface between the block and the loading plate. In this study, we adopt "velocity-weakening" 
type friction law as a dynamical friction force. (f>. This friction law states that a dynamical friction force decreases 
as the slipping velocity increases. This friction property can be observed in rock-fracture experiments jsil l and is 
formulated mathematically by Carlson et al. \\\ , namely, 

(-00,1] U = 0, 

{l + 2a[?7/(l-a)]} 



It is easy to simulate this friction law so that this formulation has been often used in previous works 
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28l |. The friction function can be characterized by two parameters, a and a. a is the difference between the 



maximum friction force {— 1) and the dynamical friction force at u = (= '/'(O)). a is represented the decrement of 
the friction force, (j). If a ~ 0, (j) is constant, 1 — <t. When a — s- oo, (f) decreases quickly to 0. v is the normalized plate 
velocity and is very small parameter. Thus, we set 1^ = when an event happens. This assumption guarantees the 
condition that no other event occurs during an ongoing event. 

In this study, we place 50 x 50 blocks on the (x, y) plane and simulate equation ^ and ([3]) under the free boundary 
condition by using the 4th-ordcr Runge-Kutta method. Initial configurations of all blocks have small irregularities. 
10^ order of events after some periods when the initial randomness effect cannot be influenced are used. We study 
the inter-occurrence time statistics by selecting lx,ly and a for tunable or control parameters, while v = 0.01 and 
cr = 0.01. This work is another version of the previous reports 



III. RESULTS AND DISCUSSION 

In this work, the slip of a block is considered as an earthquake. An earthquake occurs when a block slips for 
the first time during an event. The inter-occurrence time is defined as the time interval between successive events. 
For example, the nth inter-occurrence time can be described by t„ = t'n+i — Kn where t'„ and are the nth 
and the n + 1th earthquake occurrence time, respectively. In order to compare our results with natural earthquake 
inter-occurrence time, we introduce scaled inter-occurrence time, r' = r/f, where f is the normalized scaled time, f 
is set at 1.0 for the model analysis and at 1000 [s] for observation analysis. 
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FIG. 3; The probability distributions of the inter-occurrence time for different control parameters, lx,ly, and a as a function 
of the scaled inter-occurrence time. For (a), x (a = 2.5), o (q = 3.5), and □ (a = 4.5), while Ix = ^ and ly — \/3. For (b), 
Ix ~ ly — 1, Ix = ^ and ly — \/3, and Ix =2 and ly — 2\/2, denoted respectively, x, o, and □. b' is calculated from the slope of 
the dash line. All plots except for in the case of a = 3.5 in (a) and Ix ~ 1 and ly — \/3 in (b) are shifted vertically for clarity. 

A. Probability density distribution 



Probability density distributions of the inter-occurrence time, p(t') arc displayed in Fig. [3l For (a), the friction 
parameter a is changed from 2.5 to 4.5, whereas stiffness parameter is fixed at Ij. = 1 and ly = For (b), and 
ly are varied when a ~ 3.5. As shown these figures, probability distributions of inter-occurrence time exhibit the 
power law in short time scale region, 1 ^ r ^ 7. We calculate the power law exponent, denoted here P from the 
slope of the distribution; for example, 2.40 {a ~ 2.5), 1.96 {a ~ 3.5), 1.78 {a = 4.5) for (a). In the case of a = 3.5 
the distribution shows the power law, 1 ^ r ^ 20, so that the system exhibits a critical state approximately (type 
A). For a = 2.5, the probability of a long inter-occurrence time region, r ^ 10, is enhanced more than expected by 
the power law decay, hereafter referred to as type B. On the other hand, as when a = 4.5, the probability is less 
than predicted by the power law (type C). As for (b), /? increases when l^ and ly are enhanced, such as /3 ~ 1.94, 
1.96, and 2.71 for {l^ = ly = 1), {Ix = 1 and ly ~ \/3), and (l^ ~ 2 and ly ~ 2A/2)DThe forms and trends of the 
distributions can be categorized into three types, type A (^j, = 1 and ly ~ \/3); type B (/^r = 2 and ly = 2\/2); and 
type C (Ix = ly = 1). It is found that the form of the distribution in a long time region and the power law exponent /3 
depend on the dynamical parameters, J^, ly, and a. Our findings are different from those of previous works studying 
of the ID stick-slip model Q [22 , 23 , 24 1 . 



B. Survivor function 



Up to now, we have been discussing statistical properties of the survivor function, D{t) = prob (t > r) = 1 — 
-F(r), where F(t) is the cumulative distribution. Abe and Suzuki have analyzed the Japan and Southern California 
earthquake catalogs and found that the survivor function of the inter-occurrence time can be described by the power 



law as 
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D(t') = -, (4) 

^ ^ (l + eT)7' ^ ' 
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FIG. 4: The survivor functions of inter-occurrence time statistics. For (a), the simulation is performed 



%/3, and 



a = 3.5. For (b) the survivor distributions for the optimal case of the model, Japan earthquakes, and California earthquakes 
are shown. The fitting parameters are calculated as g = 1.06 and tq = 2.55 for the model, g = 1.05 and ro = 3.16 for Japan, 
and q= 1.13 and ro = 3.44 for California. 



where 7 and r are parameters. This distribution function can be written as 

D(t') = e,(-T/To) = [(1 + (1-(?)(-t/to))^] + , (H+ =max[0,a]), 



(5) 



where q and tq are positive constants and arc related to 7 and e: 7 = 1/('Z — 1) and e = (g — 1)/to. epfx) is the 



so-called g-exponential distribution derived from the non-additive statistical mechanics proposed by Tsallis [30|. In 
this work, we select the power law distribution defined by Eq. ([5]) for the ideal survivor distribution function, and 
then optimize the control parameters, l^^ ly and a. 

We present the survivor function of the inter-occurrence time in the case oi Ix ^ ^,ly ^ -v/S, and a ~ 3.5 in Fig. [4] 
(a). Plots and the dashed line correspond to the numerical data and the ideal curve, respectively. Varying the control 
parameters, lx,iy, and a systematically, we can find the optimal case of the control parameters: Z^. = 1, ly = \/3, and 
a ~ 3.5. For an optimal case, the fitting parameters of the survivor functions are estimated to be (7 = 1.06, tq — 2.55 
and the correlation function p yields 0.986. Figure 3] (b) shows the inter-occurrence time statistics obtained from 
the model (optimal case), Japan earthquakes, and California earthquakes. Therefore, it is concluded that the inter- 
occurrence time can be reproduced semi-qualitatively for the optimal case. 

Finally, magnitude distributions derived from the model in the case for = l,ly = a/3, and a = 3.5 
and from Japan and California earthquake catalogs are shown in Fig. O (a) and (b), respectively. Note 



that we use the JMA catalogs: 'http://kea.eri.u-tokyo.ac.jp/tseis/jmal/' for Japan and the NCEDC catalogs: 



http://www.ncedc.org/ncedc/catalog-search.htm]' for California for the period 2001/01/01-2006/08/31. It should 
be noted that the JMA catalog lacks earthquake data whose magnitude is less than 2.0. Comparing Fig. [5] (a) with (b), 
we demonstrate that the model can reproduce the GR law, which we explained previously: the power law distribution 
with exponent, 6- value 1.0. We found that in the case of optimal parameters, the inter-occurrcncc time statistics and 
the GR law can be extracted from the 2D stick-slip model [25 1. 
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IV. CONCLUSION 



In this study we numerically investigated the statistics of earthquake interval times, the inter-occurrence time based 
on the 2D spring-block model. Inter-occurrence times arc interval times between earthquakes on all faults in a region. 
It is found that intcr-occurrcncc time statistics depend on the control parameters, l^, ly, and a charactering the model. 
Ix and ly are stiffness parameters, whereas a expresses the decrement of dynamical friction. The probability density 
distributions of the inter-occurrence time show the power law in the short-time region. For the long time region, the 
distributions could be classified into three types: power law behavior (type A), broad peak structures (type B), and 
exponential cutoffs (type C). Then, we restricted the control parameters so as to reproduce the inter-occurrence time 
statistics in nature; the survival function of the inter-occurrence time revealed the (/-exponential distribution with 
q > 1. The optimal parameters are estimate to be /j, = 1, ly = \/3, and a = 3.5. 

In the case of the optimal parameters, the magnitude distribution shows the power law with exponent 1.0. This 
power law distribution is similar to the GR law, and the exponent corresponds to the 6-value, which is characterized 
by the GR law. Hence, we demonstrate that in the optimal case, Ix = l,ly = VS, and a ~ 3.5, the model can 
reproduce the GR law and the inter-occurrence time statistics both simultaneously and spontaneously. 

We acknowledge that the stick-slip (block-slider) model is highly simplified so that many effects playing important 
roles in fault dynamics have been neglected. However, it is shown that this model is useful for the extraction of the 
statistical properties of earthquakes because the inter-occurrence time statistics and the GR law can be extracted. 
This work is a first step toward studying the origin of the statistics of time intervals. We hope to extend our work by 
focusing on recurrence time statistics and comparing them with natural recurrence time statistics. 
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FIG. 5: The magnitude distribution obtained from the model in (a) and from real earthquakes in (b). The b- value of the GR 
law is calculated from the slope of the distribution: 1.05 (model), 0.83 (Japan), and 1.02 (California). 
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